{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "1f925241-060a-45d3-85d6-d5f41ecb4cc2",
   "metadata": {},
   "source": [
    "# Validaton of the spherical quadrature equations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f569718c-a9bc-47e5-80fb-255f70044c7b",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy\n",
    "import numpy.matlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "890da482-b7fc-4ee6-aab0-910b9936732f",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rquad(N, k):\n",
    "    k1 = k+1; k2 = k+2; n = numpy.arange(1, N+1); nnk = 2*n+k\n",
    "    A = numpy.insert(numpy.matlib.repeat(k**2, N) / (nnk*(nnk+2)),0,k/k2 )\n",
    "    n = numpy.arange(2, N+1); nnk = nnk[1:N+1]\n",
    "    B1=4*k1/(k2*k2*(k+3)); nk=n+k; nnk2=nnk*nnk\n",
    "    B=4*(n*nk)**2/(nnk2*nnk2-nnk2)\n",
    "    ab = numpy.column_stack((A,numpy.concatenate(([(2**k1)/k1], [B1], B))))\n",
    "    s = numpy.sqrt(ab[1:N,1])\n",
    "    [X,V] = numpy.linalg.eig(numpy.diag(ab[0:N,0],0)+numpy.diag(s,-1)+numpy.diag(s,1))\n",
    "    I = numpy.argsort(X)\n",
    "    X=numpy.sort(X)\n",
    "    V = V[:,I]\n",
    "    x=(X+1)/2\n",
    "    w=(1/2)**(k1)*ab[0,1]*V[0]**2\n",
    "    return x, w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "a6c238b7-5f86-47fc-bf90-c3f59a66ce24",
   "metadata": {},
   "outputs": [],
   "source": [
    "def spherequad(nr,nt,np,rad):\n",
    "    r, wr = rquad(nr,2);         # radial weights and nodes (mapped Jacobi)\n",
    "\n",
    "    if rad == float('inf'):        # infinite radius sphere\n",
    "        wr=wr/(1-r)**4           # singular map of sphere radius\n",
    "        r=r/(1-r) \n",
    "    else:                        # finite radius sphere\n",
    "        wr=wr*rad**3             # Scale sphere radius\n",
    "        r=r*rad\n",
    "    x,wt =rquad(nt,0)\n",
    "    t=numpy.arccos((2*x-1)); wt=2*wt       # theta weights and nodes (mapped Legendre)\n",
    "    p=2*numpy.pi*numpy.linspace(0, np-1, np)/np         # phi nodes (Gauss-Fourier)\n",
    "    wp=2*numpy.pi*numpy.ones(np)/np        # phi weights\n",
    "    rr,tt,pp = numpy.meshgrid(r,t,p)   # Compute the product grid\n",
    "    \n",
    "    r = rr.flatten('F')\n",
    "    t = tt.flatten('F')\n",
    "    p = pp.flatten('F')\n",
    "\n",
    "    wt = wt[:, numpy.newaxis]\n",
    "    wr = wr[:, numpy.newaxis]\n",
    "    wp = wp[:, numpy.newaxis]\n",
    "    w=numpy.reshape(\n",
    "        numpy.dot(numpy.reshape(numpy.dot(wt,wr.transpose()),(nr*nt,1),'F'),\n",
    "        wp.transpose()),(nr*nt*np,1),'F')  \n",
    "    w = w.reshape(-1)\n",
    "    return r, t, p, w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b64f4920-9abf-442c-a756-65a110ecc2b7",
   "metadata": {},
   "outputs": [],
   "source": [
    "r,t,p,w = spherequad(4,5,6,2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6e8510ab-c0ca-49f2-a8d3-0aec5f782710",
   "metadata": {},
   "outputs": [],
   "source": [
    "t"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "25384659-7eec-4cd7-9565-75d305d78590",
   "metadata": {},
   "outputs": [],
   "source": [
    "p"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0bc26ed3-83e0-4595-9f2e-d455631d168e",
   "metadata": {},
   "outputs": [],
   "source": [
    "w"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "43111db4-ecfb-4318-97fd-14b10c237db2",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "270fae73-4468-47da-a38b-d0e9cc747c07",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}
